#Figure A10

#directory
#set working directory to the path PvP_Replication
#setwd("~/PvP_Replication")

#libraries
library(sf) #spatial data
library(tidyverse) # data manipulation
library(ggtext) #text labels

#read in basemap of Colombia for use with all four maps
col_shapefile <- st_read("PvP_data/country_shapefile/ne_10m_admin_0_countries.shp") %>% 
  dplyr::select(NAME, geometry) %>%
  filter(NAME == "Colombia")

#step 1 plot (top left pane)

#read in coca cluster data
coca_clusters <- st_read("PvP_data/coca_clusters/coca_clusters.shp")
coca_clusters <- coca_clusters %>% mutate(clusterfactor = as.factor(cluster))

#r crop the two shapefiles to allow for plotting; full file too big to plot
coca_clusters_crop <- st_crop(coca_clusters, xmin = -77.79753785173181768, ymin = 7.32443275890097745, 
                                  xmax = -76.20867611536959885, ymax = 9.05506161988810909)
col_shapefile_crop <- st_crop(col_shapefile, xmin = -77.79753785173181768, ymin = 7.32443275890097745, 
                            xmax = -76.20867611536959885, ymax = 9.05506161988810909)

#plot the cropped shapefiles, coloring by the cluster variable
coca_cluster_plot <- ggplot() +
  geom_sf(data = col_shapefile_crop, aes(fill = "Colombia"), color = NA) +
  geom_sf(data = coca_clusters_crop, aes(color = clusterfactor)) +
  scale_fill_manual(values = "#E0E0E0", name = "Colombia") +
  coord_sf(expand = FALSE) +
  guides(colour=FALSE, fill = FALSE) +
  theme_void()

#save step 1 plot

ggsave(plot = coca_cluster_plot, filename="PvP_plots/coca_cluster_plot.pdf", width = 8.5,
       height = 11, units = "in")


#step 2 plot (top right pane)

#read in ports and border crossing dataset
exit_points <- read.csv("PvP_data/exit_points_18.csv")

#convert to point geometry
exit_points_sf <- st_as_sf(exit_points, coords = c("to_long", "to_lat"), crs = "+init=epsg:4326")

#plot exit points, coloring by type  
exit_points_plot <- ggplot() +
  geom_sf(data = col_shapefile, fill = "#E0E0E0", color = NA) +
  geom_sf(data = exit_points_sf, aes(color = end_type), size = 1.5) +
  scale_color_manual(values = c("port" = "darkblue", "border" = "darkred"), labels = c("Border Crossings", "Ports")) +  
  labs(title = "",
       color = "") +
  theme_void() +
  theme(legend.position = "bottom")

#save step 2 plot
ggsave(plot = exit_points_plot, filename="PvP_plots/exit_points_plot.pdf", width = 8.5,
       height = 11, units = "in")

#step 3 plot (bottom left pane)

#read in simulated trafficking routes
trafficking_routes <- st_read("PvP_data/trafficking_routes/trafficking_routes.shp")

#plot simulated trafficking routes
trafficking_route_plot <- ggplot() +
  geom_sf(data = col_shapefile, fill = "#E0E0E0", color = NA) +
  geom_sf(data = trafficking_routes, color = "darkred", lwd = 0.25) + 
  theme_void()

#save result
ggsave(plot = trafficking_route_plot, filename="PvP_plots/trafficking_route_plot.pdf", width = 8.5,
       height = 11, units = "in")


#step 4 plot (bottom right pane)

#select simulated routes to plot
sample_routes <- trafficking_routes %>% filter(start_id %in% c(114, 100) & end_id == 8)

#select exit point to plot
sample_exit <- exit_points_sf %>% filter(end_id == 8)

#extract intersections of routes to identify simulated hub for plot
sample_hub <- sample_routes %>% 
  st_union(by_feature = FALSE) %>% 
  st_line_merge() %>% 
  st_cast("LINESTRING") %>% 
  st_intersection() %>% 
  st_as_sf() %>% 
  mutate(geomtype = st_geometry_type(x)) %>% 
  filter(geomtype == "POINT") %>% 
  st_union(by_feature = FALSE) %>% 
  st_cast("POINT") %>% 
  st_as_sf()

#add label for plot
sample_hub$name <- "Hub"

#plot all five layers
sample_hub_plot <- ggplot() +
  geom_sf(data = col_shapefile_crop, aes(fill = "Colombia"), color = NA) +
  geom_sf(data = coca_clusters_crop, aes(color = clusterfactor)) +
  geom_sf(data = sample_routes, color = "darkred") +
  geom_sf(data = sample_exit, size = 3, color = "darkblue") +
  geom_sf(data = sample_hub, size = 2.5, color = "black", shape = 15) +
  geom_richtext(data = sample_hub,
                aes(label = name, geometry = x), 
                stat = "sf_coordinates",
                size = 9 / .pt,
                colour = "black",
                fontface = "bold",
                label.colour = "black",
                fill = "#E0E0E0",
                label.r = unit(0.15, "lines"), nudge_x = -0.1
  ) +
  scale_fill_manual(values = "#E0E0E0", name = "Colombia") +
  coord_sf(expand = FALSE) +
  guides(colour=FALSE, fill = FALSE) +
  theme_void()

#save result
ggsave(plot = sample_hub_plot, filename="PvP_plots/sample_hub_plot.pdf", width = 8.5,
       height = 11, units = "in")

